function [dof_map, V, T, u1, u2, err] = ex3(n_round, FE_Type, FE_Order, Young, nu, method)
% function [dof_map, V, T, u1, u2, err] = ex3(n_round, FE_Type, FE_Order, Young, nu, method)
%
%  Solve planar elasitic problem on domain L: [-1,1]x[-1,1]/[0,1]x[0,1] 
%  with analytic solution given by fun_u
%  the right hans side is then given by fun_f
%
%  see the example in paper of M. Ainsworth and B. Senior:
%  Comput. Methods Appl. Mech. Engrg. 150(1997), 65-87.
%
%  This is also for validation!!
%
%  Test pass on Aug. 11th, 2012
%
%  Author: Dr. Xian-Liang Hu
%

if nargin < 1
    FE_Type = 'BB';  % default FE basis.
    FE_Order = 2;   % default polynomial order.
end

if nargin < 3
    % default physical parameters for elastic problems
    Young = 1; 
    nu = 0.3;
end

if nargin < 5
    method = 2;  % using FEM style by default
end

% %% for output dofs v.s. errors
% dof_count = zeros(10,1);
% err_L2 = zeros(10,1);
% err_inf = zeros(10,1);

Quad_Order = 13;

% generate initial mesh
[p, e, t] = initmesh('lshapeg','Hmax',0.35,'Hgrad',1.99);
V = [p(1,:)' -p(2,:)']; T = t(1:3,:)';
[T, E, ET, TE] = build_fem_mesh(V, T);
[V, T, TE, E, ET] = refine_L_shape(V, T, TE, E, ET, 3);
V = V*sqrt(2);  % scale
x = V(:,1); y = V(:,2); theta = 45/180*pi;
V(:,1) =  cos(theta)*x + sin(theta)*y;
V(:,2) = -sin(theta)*x + cos(theta)*y;
for round = 1:n_round
    [V,T, TE, E, ET] = refine_mesh_uniform(V, T, TE, E, ET);
end

fprintf('There are %d triangles.\n', size(T,1));


% [dof_map, n_dof] = distribute_dof(T, TE, ET, FE_Order, 0);
% [dof_x, dof_y] = build_dof_coordinate(dof_map, V, T, FE_Order);
%   
% [K11, K12, K22, b1, b2] = elastic_T3(V, T, Quad_Order, FE_Type, FE_Order, @fun_f, Young, nu);
% [A, b] = assemb_system(dof_map, K11, K12, K22, b1, b2); 
% 
% % find the sigularity:
% dof_singular = find(abs(dof_x) < 1e-6 & abs(dof_y) < 1e-6);
% 
% % mode 1: y-direction displacement = 0;
% dof_singular = dof_singular + n_dof;
% A(dof_singular,:) = 0;
% b(dof_singular) = 10;
% A(:,dof_singular) = 0;
% A(dof_singular, dof_singular) = 1;
% % 
% % % mode 2: x-direction displacement = 0;
% % A(dof_singular,:) = 0;
% % b(dof_singular) = 0;
% % A(:,dof_singular) = 0;
% % A(dof_singular, dof_singular) = 1;
% 
% c = A\b;
% u1 = c(dof_map);
% u2 = c(dof_map + n_dof);


%%%%%%%%
%  
bdr_Dirichlet = find(ET(:,2)==0);
bdr_Neumann = [];

[dof_map, u1, u2, t] = elastic_fem(V, T, TE, ET, FE_Type, FE_Order, Quad_Order, method, ...
                                     Young, nu, @fun_f, bdr_Dirichlet, bdr_Neumann, @fun_u);


% finite element solutions:
disp_d = 1; tri_temp = template_mesh_tri(disp_d);
[u_h, Tris1, Points1] = fe_solution_bb(V,T, u1, FE_Order, tri_temp, disp_d);
[v_h, Tris2, Points2] = fe_solution_bb(V,T, u2, FE_Order, tri_temp, disp_d);

% analytic solutions
[uu, vv] = fun_u(Points1(:,1),Points1(:,2), Young, nu);

% cal the L_inf error
erru = max(max(abs(u_h - uu)));
errv = max(max(abs(v_h - vv)));

% return the maximum
err = max(erru, errv);

% % visualization and print information
subplot(2,2,1); trisurf(Tris1, Points1(:,1), Points1(:,2), u_h - uu);  %plot error for u
subplot(2,2,2); trisurf(Tris1, Points1(:,1), Points1(:,2), v_h - vv);  %plot error for v
subplot(2,2,3); fac = 0.1; ZERO=zeros(size(Points1,1),1);
trisurf(Tris1, Points1(:,1)+u_h*fac, Points1(:,2)+v_h*fac,ZERO);
view(0,90);axis tight; axis equal;
subplot(2,2,4); fac = 0.02; ZERO=zeros(size(Points1,1),1);
trisurf(Tris1, Points1(:,1)+u_h*fac, Points1(:,2)+v_h*fac,ZERO);
view(0,90);axis tight; axis equal;
fprintf('The inf norm or error(u) = %e,  error(v) = %e.\n',  erru, errv);

end

% %%%%%%%%%%%%%%%%%%%%%%%%
% % mode 1 boundary condition u_2(0,0) = 0;
% function [constrains, u_1, u_2] = mode_1_boundary()
% 
% end
% 
% %%%%%%%%%%%%%%%%%%%%%%%%
% % mode 1 boundary condition u_1(0,0) = 0;
% function [constrains, u_1, u_2] = mode_2_boundary()
% 
% end


%%%%%%%%%%%%%%%%%%%%%%%%
%  analytic solution, also used as boundary function
% 
function [u_1, u_2] = fun_u(xx, yy, E, nu, varargin)
    kappa = 3 - 4*nu;
    G = E/(2 + 2*nu);
    r = sqrt(xx.*xx + yy.*yy);
    idx = (r < 1e-10);
    theta = atan(yy./(xx+1e-14));

    % for mode 1:
     lambda = 0.5444837367825;  Q = 0.5430755788367;
     a = kappa - Q*(lambda + 1); 
     b = kappa + Q*(lambda + 1); 
     u_1 = 0.5/G*r*lambda.*(a*cos(lambda*theta)  - lambda*cos(lambda - 2)*theta);
     u_2 = 0.5/G*r*lambda.*(b*sin(lambda*theta) + lambda*sin(lambda - 2)*theta);
    
    
%     % for mode 2:
%      lambda = 0.9085291898461;  Q = -0.2189232362488;
%      a = kappa - Q*(lambda + 1); 
%      b = kappa + Q*(lambda + 1); 
%      u_1 = 0.5/G*r*lambda.*(a*sin(lambda*theta)  - lambda*sin(lambda - 2)*theta);
%      u_2 = -0.5/G*r*lambda.*(b*cos(lambda*theta) + lambda*cos(lambda - 2)*theta);

    
    u_1(idx) = 0; u_2(idx) = 0;  % displace at the original is 0
end


%%%%%%%%%%%%%%%%%%%%%%%%
%  The right hand side function
% 
function [f_1, f_2] = fun_f(xx, yy, E, nu, varargin)
    f_1 = zeros(size(xx));
    f_2 = f_1;
    
end